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The Food and Drug Administration (FDA) 
is conducting research on developing 
reference lung cancer lesions, called 
phantoms, to test computed tomography 
(CT) scanners and their software. FDA 
loaned two semi- spherical phantoms 
to the National Institute of Standards 
and Technology (NIST), called Green and 
Pink, and asked to have the phantoms' 
volumes estimated. This report describes 
in detail both the metrology and computa- 
tional methods used to estimate the 
phantoms' volumes. Three sets of coordi- 
nate measuring machine (CMM) measured 
data were produced. One set of data 
involved reference surface data measure- 
ments of a known calibrated metal sphere. 
The other two sets were measurements of 
the two FDA phantoms at two densities, 
called the coarse set and the dense set. 
Two computational approaches were 
applied to the data. In the first approach 
spherical models were fit to the calibrated 
sphere data and to the phantom data. The 
second approach was to model the 
data points on the boundaries of the 
spheres with surface B-splines and then 
use the Divergence Theorem to estimate 
the volumes. Fitting a B-spline model to 
the calibrated sphere data was done 
as a reference check on the algorithm 
performance. It gave assurance that the 
volumes estimated for the phantoms would 



be meaningful. The results for the coarse 
and dense data sets tended to predict the 
volumes as expected and the results did 
show that the Green phantom was very 
near spherical. This was confirmed by 
both computational methods. The spherical 
model did not fit the Pink phantom as well 
and the B-spline approach provided a 
better estimate of the volume in that case. 
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1. Introduction 

Tomography is a method of imaging a single slice of 
the body. Modern computed tomography (CT) is a 
medical imaging method that uses tomography, but also 
employs digital image processing techniques, to gener- 
ate three dimensional images built from a large 
sequence of two-dimensional x-ray images made 
around a single axis. CT has shown promising results in 
detecting lung cancers at more operable stages, when 
survival is better [1]. 

The Food and Drug Administration (FDA) is con- 
ducting research on developing reference cancer 



lesions, called phantoms, to test CTs and their software. 
Two samples were loaned to NIST to estimate volumes. 
The material of the phantoms simulates lung cancer 
material. The phantoms can be inserted into a simulat- 
ed body torso for CT scans. The two phantoms are 
shown in Fig. 1 . Although they seem spherical, they are 
slightly non- spherical. The larger one on the right is 
referred to as the Green phantom and the one on the left 
is called the Pink phantom due to the material colors. 

One experimental approach to estimate the volume of 
the phantoms would be to use an Archimedes test in 
which the phantoms would be immersed in a liquid bath 
in a well calibrated container with fine measurement 



149 



Volume 115, Number 3, May- June 2010 

Journal of Research of the National Institute of Standards and Technology 




Fig. 1. Two simulated lung cancers, called phantoms. 

gradations to determine liquid displacement. However, 
in the case of these phantoms, the material used to 
manufacture these phantoms was porous and the phan- 
toms would have to be coated. This, of course, would 
affect the "ground truth" volume estimate. As a result, 
the approach chosen for this study was based on a fun- 
damental theorem in calculus, called the Divergence 
Theorem (see Taylor [2]), an analogue of Green's 
Theorem in two dimensional space. In the Divergence 
Theorem a volume integral is shown to be equal to a 
surface integral. Therefore, we surmised that, if a 
model of the surfaces of the phantoms could be devel- 
oped, the Divergence Theorem would help to estimate 
their volumes. In order to develop a surface model we 
needed to obtain data about the surfaces of the phan- 
toms. This was done using a coordinate measuring 
machine (CMM) in the Manufacturing Engineering 
Laboratory (MEL) at NIST. This machine produced a 
set of (x, y, z) points on the surface of each phantom. 
The data were then transformed to spherical coordi- 
nates and modeled using a set of basis functions, called 
B-splines. After fitting, the B-spline model was 
employed to generate a grid of values on the surface. 
These values were used to form surface triangles that 
were then used to compute the necessary surface 
integrals and finally the volume. The quality of the 
volume estimates depended on the surface grid sizes, as 
will be clear from the discussion below. The method of 
B-spline surface modeling is not new, in that it was 
suggested in the book by Dierckx [3]. However the 
application in the current case, that uses the Divergence 
Theorem, appears to be new. A reader can also consult 
the Dierckx references [4] and [5]. 

The paper is divided as follows. Section 2 describes 
how a volume can be computed using the Divergence 



Theorem. Section 3 describes the surface point genera- 
tion experiments using the CMM. Section 4 introduces 
the two methods of data modeling used to estimate the 
phantoms' volumes. The first uses a spherical model. 
The second uses the B-splines as a basis for a least 
squares fit with regularization to model the phantoms' 
surface data. The fitted functions were then used to 
generate the grid data necessary to apply the 
Divergence Theorem. Section 5 presents the specific 
surface triangulation used to implement the Divergence 
Theorem. Section 6 describes the volumetric results 
from applying the spherical and B-spline models to the 
CMM data. A summary is given in Sec. 7. 



2. Volume Estimation by the Divergence 
Theorem 

In this section we state the Divergence Theorem and 
indicate how it can be used to estimate the volume of a 
polyhedron. This provides the motivation for the need 
to determine surface data points from the phantoms and 
then to build a surface model that is used to create data 
points at surface triangulation vertices. As the number 
of triangles was increased, the volume estimates were 
expected to tend to a fixed value. As will be seen, this 
expectation is confirmed below. 

2.1 Divergence Theorem in 3-D and Volume 
Computation 

Let Q be a simply connected region in three dimen- 
sional space and T the surface boundary. Then 



[[[ V • F dxdydz = [ [ F -h do , 



(1) 



Q. 



where F (x, y, z) = (F^ (x, y, z), F2 (x, y, z), F^ (x, y, z)) 
is a differentiable vector field. 



VF = 



dF dF dF 



+ 



2^-3 



(2) 



dx dy dz 



is the divergence of the vector field, fi is an outward- 
pointing unit normal vector to F, and da is an infinites- 
imal element of the surface. 

If we select F (x, y, z) = (l/3)(x, y z), then y'F=l 
and we can write 



dxdydz = —\\n(x,y,z) da 



(3) 



Q 
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Note that the left hand side is simply the volume of the 
region Q. The 1/3 factor comes from the definition of 
F (x, y, z) so that V- F = 1 . Now, if T is approximated 
by disjoint polyhedra (planar surface patches), F,, then 



A^ 



r = Ur, 



(4) 



i=\ 



and 



1 1 

-JJ^(x,j;,z)(icT=-^JJ/2. -{x,y,z) da, (5) 



where n^ is the unit outer normal to F,. Here we will 
model the surface patches by planar facets and, in 
particular, triangular facets. The plane for F, is given by 
rii • (x, y, z) = c,, where c^ is a constant associated with 
each triangular facet. The sum of the integrals over the 
facets thus reduces to 

1^ Cc 1^ rr 1^ 

-Z J J n,ix,y,z)da= -2^ J J q ^^=tZ ^r Area(F.). 

-^ i-l r 2-1 r 2=1 

i i 

(6) 

This method of computing the volume of an object 
from a surface integral can be found in Schneider and 
Eberly [6]. 

2.2 Area of a Planar Polyhedron in 3-D 

In order to estimate the surface area of a facet, as 
used in (6), we will describe the process of computing 
the area of a planar polyhedron in space by use of 
Stokes' Theorem and then we will particularize it to a 
planar triangle in space. Here the assumption is made 
that there is surface data available at the planar triangle 
vertices. Stokes' Theorem states that if C is a piecewise 
smooth boundary curve, oriented positively, of a sur- 
face F, and if F is a differentiable vector field defined 
on F and /? is a unit normal satisfying the right-hand 
rule relative to the boundary orientation, then 



^^{VxF)nda=^FdR, 



(7) 



c 



where the curl oiF{x, y, z) = (Fi(x, y z), F2(x, y, z), 
F3 (x, y z)) is 



VxF = 



/ 



V 



3^3 dF^ dF, 3^3 dF^ dF, 



dy dz dz dx dx dy 



\ 



J 



(8) 



and dR -(dx, dy, dz). If we set F{x, y, z) =(1/2) 
n X (x, y, z) then V xF = h and ri-h=lso that 

Area (F) = [ [ da = — (p (/z x(x, y, z)) • ( dx, dy, dz) , 



c 



= — d) /z • ((x, y, z) x(dx, dy, dz)) , 



c 



-jn-iixit),yit),zit))x 



c 



(x'(t),y'it),z'it)))dt, 
^j>n-iyit)z'it)-zit)yXt), 



(9) 



C 



z(t)xXt)-x(t)zXt), 
x{t)y\t)-y{t)xXt))dt. 

The factor 1/2 comes from the definition of F{x,y,z) 
so that (VxF) = h. 



We will now specialize the Stokes formula to find the 
area of a spatial triangle, y, in terms of its three vertices 
oriented positively. Let the three vertices, in positive 
orientation, be specified by v^ = (x^, y^, z^), 
V2 = (x2, y2, ^2)5 ^3 ^ fe? J^35 ^3)- Parameterize the bound- 
ary of the triangle as follows. For /g[0, 1], 
let v = Vi + /(v2-Vi) = (xi + /(x2-Xi), yi + t(y2-yi), 
Zi + /(Z2-Z1)). For /g[1,2] let v = V2 + (/-l) 

(V3 - V2) = (x2 + (/ - 1 )(x3 - X2), y2 + (t-i )(y3 - yi), 

Z2 + {t - l)(z3 - Z2)). Finally, for / g [2, 3], put 
V = v^ + {t - 2)(vi - V3) = (X3 + {t - 2)(xi - X3), 
y3 + {t-2){y^-y^), Z3 + (^-2)(zi-Z3)). The area of 
the spatial triangle, in terms of the parameterized 
boundary, can be written as 



Area(y) =^^«-(y(Oz'(0 -z(0/(0, 



C 



z{t)x{t)-x{t)z\t), 
x{t)y\t)-y{t)x\t))dt, 

^h-(^\l{y{t)zXt)-z{t)yXt)) dt, 

\ (z(t)xXt)-x(t)zXt)) dt, 

J 

fjxit)y'it)-yit)xXt))dt) 



(10) 
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By straightforward integration over each of the 
parameterized segments it is easy to show that 

{x,y^-x^y,)+{x^y^ -x^y^)^x^y, -\y^)) . 

(11) 

The normal vector to the oriented triangle can be com- 
puted as follows. Let v =(^2 -x^ ^2 ~ Ji' ^2 ~^i) ^^^ 
w= (x^-x^, y^- yi, -^3 - -^1 ) • Then n = vxw = ((^2 ~ 

yi)(^3 -^i)-(^2 -^1X3^3 -Ji)' (^2 -^X^ -^) -(^ - 

^1X^3 -^iX (^2-^X^3-3^1) -(3^2 -JiX^ -^)).If we 
set n, =(y^-y,)(z^-z^)-(z^ -z^){y^-y^\n^ =(^2- 

zj(x3 - X,) - (X2 - ^)(Z3 -z,\ n^ = (x, - j^)(y^ - 



yi) - te - 7iX^ - ^X then 
and n = n/ 



n 



2 , 2 , 2 
n^ + ^2+ ^3 



n 



3. Surface Metrology of the Artifacts 

In this section we discuss the experimental method 
used to obtain surface data for the three objects used in 
this study. As a test object for the modeling process, a 
well calibrated sphere was selected. This object, along 
with the FDA phantoms formed the study artifact set. 
Data points, in the form of {x, y z) coordinates, were 
created by the probing action of the CMM. Figure 2 
shows the CMM system used to measure the artifacts. 
The system is computer controlled and touches an 
object to be measured at programmed points in order to 
produce {x, y z) values at the probed points. 

Figure 3 shows one of the phantoms, held by a 
device called a vacuum chuck, as it is being touched by 
the CMM probe. The probe itself can be programmed 
to approach an object at various angles. In the back- 
ground of the figure is a high quality reference steel 
sphere. Before an object is probed the reference sphere 
is measured to determine the effective diameter of the 
CMM probe tip. This reference sphere is separate from 
the calibrated sphere used to test the modeling process. 
The difference between the known diameter of the ref- 
erence steel sphere and the apparent diameter of the 
measured reference metal sphere gave the calibrated 
effective probe diameter. 




Fig. 2. Computer controlled CMM. 




Fig. 3. Lung nodule phantom on the vacuum chuck. 
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The FDA phantoms were created by a molding 
process and the mold marks on both spheres were 
visible. These mold marks were used to align the 
phantoms with the coordinate system of the CMM. The 
marks were laid out like lines of latitude, leading to the 
use of a natural nomenclature of latitude and longitude 
like those of the Earth. For each phantom, the North Pole 
was chosen to be the one with darker, deeper, or more 
obvious mold marks. 

For the purpose of understanding the measurement 
process we will describe the physical fixture positioning 
of the phantoms. They were held by a vacuum chuck (see 
Fig. 3) to minimize distortion and reduce the chance of 
damaging the spheres' surfaces. As Fig. 3 shows, the 
chuck has a shallow cone to hold the phantoms. The 
phantoms contacted the cone around a circle of latitude 
at about - 45°. This was measured from the equatorial 
circle around the middle of the phantoms. For example, 
the point at the top of the phantoms would be at + 90°. 
The wall vacuum of the cone was strong enough to hold 
the spheres sufficiently that they did not move signifi- 
cantly, as shown by the repeatable results from run to run 
at the 1 |im level. The setup of the phantoms in the 
vacuum chuck was accomplished by eye alignment 
using the visible mold marks and minor imperfections as 
guides to the eye. The expanded uncertainty of alignment 
was estimated to be approximately ± 2° for the vertical 
angle (i.e., keeping the equator horizontal) and ± 4° for 
the azimuthal angle. For a presentation of the guidelines 
for how expanded uncertainty of parameters is comput- 
ed see Taylor and Kuyatt [7]. Essentially, the guideline 
for estimating the expanded uncertainty involves apply- 
ing a multiplier, k=2, times the square root of the vari- 
ance about a predicted parameter value (i.e., two times 
the standard error). For a full discussion of confidence 
intervals for parameter estimation see Draper and 
Smith [8]. 

Three sets of measurements were planned for each of 
the phantoms. In each, points were probed on only a 
hemisphere at a time and later the data sets were post- 
processed to form a spherical data set. In the first set of 
measurements the North Pole was set as the top point. 
A set of points was programmed to probe the top 
hemisphere down to the equator. The phantom was then 
re-positioned in the chuck so that the South pole was 
the top point. The same hemispherical points were 
probed. The rotation was done in such a manner that the 
mold marks representing the longitudinal lines were 
kept as aligned as possible. Post-processing of the data 
associated the correct signs with the measured coordi- 
nates relative to the CMM coordinate system. The third 
measurement involved re-positioning the phantoms so 



that the equatorial circle was vertical and the North- 
South axis was horizontal. Again two hemisphere sets 
of probe points were measured. This position was not 
feasible for the Pink phantom in one of the experiments 
described below. 

Visually, the pink sphere was noticeably out-of- 
round, in the shape of an oblate spheroid. A dial caliper 
gave diameter measurements given in Table 1. The 
uncertainty of caliper measurements on hard steel 
surfaces is about 0.1 mm, and is estimated to be about 
0.3 mm on the sample spheres due to the potential that 
the contact force would distort the soft surfaces of the 
spheres. All estimated uncertainties were k=2 expanded 
uncertainties. 

Table 1. Dial Caliper Measurements of Phantoms' Diameters 



Equatorial 



Polar 



Pink 
Green 



20.0 mm 
20.0 mm 



18.4 mm 
20.1 mm 



Two probing experiments were performed on each of 
the phantoms and the calibrated sphere. They created 
what we will call a coarse data set and a dense data set. 
For the coarse data set the plan was to measure each 
phantom on the CMM three times in each position, 
with 61 coordinate points per hemisphere. For the green 
sphere, each measurement set consisted of three 
separate sets of points: North pole up. South pole up, 
and prime meridian/equator intersection up. The pink 
sphere could not be held sideways in the first experi- 
ment, as the out-of-roundness prevented an effective 
vacuum seal. Therefore, a measurement data set for this 
sphere had only the North Pole up and South Pole up 
data. The third data set in this case was a re-measure of 
the North Pole up position. Each data set consisted of 
122 points. The plots in Fig. 4 and 5 show the radial 
deviation from a best- fit sphere for the full data sets. 
The figures show the radial residuals obtained by fitting 
sphere models to the Green and Pink data with an 
algorithm ordinarily used during sphere calibration 
work. In particular, they represent the residual errors 
between the distance from the fitted sphere center to the 
probed points and the fitted radius of the sphere model. 
The residuals, in the case of the Green phantom, range 
from approximately -0.1mm to +0.1 mm, whereas 
the residuals, in the case of the Pink phantom, range 
from approximately - 0.57 mm to + 0.23 mm. This 
suggests the slight non-spherical nature of the Pink 
phantom. Figure 6 shows the typical distribution of the 
probe points on a sphere. The plot is a transparency so 
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n 



CMM Pia/be Points 




Fig. 6. Sample Distribution of the CMM Probe Points on a Sphere. 



that probe points on the opposite side of the sphere are 
visible. The calibrated sphere on the shaft, with a dia- 
meter of 19.05 mm, was also measured with 122 points, 
repeated five times. It was measured in one position, 
since it was permanently mounted on a support shaft. 

For the dense data set 181 points were taken on the 
phantoms and the calibrated sphere, with five repeats in 
each position. The positions were taken the same as 
those in the first experiment. That is, the alignment was 
selected with North pole up (position 1), South pole up 
(position 2), and prime meridian/equator intersection 
up (position 3). In this case it was possible to hold the 
Pink phantom in the sideways position 3. The calibrat- 
ed sphere was also measured at 181 points with five 
repeats. 



4.1 A Spherical Model 

In this section we consider how close the data could 
be modeled by a spherical model for the data. The 
calibrated sphere, of course, could be modeled via a 
spherical model. 

In particular, let c = (c^, C2, C3, C4) and set 



F(x,y,z,c)=(x-c,y+(y-c^f+(z-c,f-c,\ (12) 

The unknown parameters c^, C2, c^ represent the 
center of the sphere and C4 is the radius. All of the data 
were measured in millimeters so that the parameters 
naturally have millimeter units. Define 



4. Modeling Methodologies 

In this section two forms of data modeling will be 
discussed. Since the phantoms seemed to be nearly 
spherical the natural tendency was first to consider 
fitting a spherical model to each phantom and estimat- 
ing the volume of the fitted spheres. However, in order 
to develop a potentially more accurate volume estima- 
tion model, the surface data was also fit using tensor 
products of B-splines and the volumes estimated by the 
Divergence Theorem. 



Oic) = ^^{(x, -cy +iy,-c,y +(z, -c,f-c,'], (13) 



i=l 



where O, the objective function, is a measure of the 
residual for the fitted sphere defined by the coefficients 
c. Since the function O is a nonlinear implicit function 
of the parameters we needed to use a nonlinear mini- 
mization algorithm to find the best fit, i.e., to solve the 
problem 



min 0(c) . 



(14) 
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There are a number of algorithms for fitting least 
squares models to data on geometric shapes (see 
Shakarji [9]). There are also various algorithms for 
minimizing general nonlinear functions, such as (13). 
All of the algorithms involve iterative minimization of 
some form. Many require computing derivatives of the 
objective function in order to generate search directions 
along which to identify a minimum. Others do not 
involve derivatives but may be somewhat slower in the 
minimization search. The algorithm selected here, 
because of the relatively few parameters involved, and 
the fact that derivatives are not required, is a form of 
polyhedron search method called the Nelder-Mead 
method (see Sauer [10]). The Newton method, requir- 
ing derivatives, was initially used to estimate the 
parameters, but the Nelder-Mead tended to produce the 
smallest value to (13). 

The Nelder-Mead algorithm works iteratively 
through steps that involve reflections, vertex exten- 
sions, and multidimensional polyhedron contractions 
until the volume of the polyhedron becomes less than a 
prescribed tolerance. The polyhedron vertices then 
enclose the minimum. The median value of the poly- 
hedron vertex values is taken as the minimum value. 

The uncertainties of the estimated center and radius 
were computed using the methods proposed in Draper 
and Smith [8] for nonlinear regression. In particular, if 
c = (ci, ^2, C3, C4), then define the nx4 matrix with 
elements 






(x,v.,z.,c), i=l,---n, 7=1,2,3,4 



(15) 



The i-th row is given by 



^/=[-2(^.-q) -2(7,-^2) -2(z.-q) -2c,]. 



(16) 



The standard error of q, s.Q.(q), is given by the 
square root of the i-th diagonal element of (Z' ZY^ s^ 
where s = 0{c) I (n-A). The expanded uncertainty is 
given by 2s.e.(c^) as defined in Taylor and Kuyatt [20]. 
The uncertainty limits of q are c, ± 2s.e. (c,). All units 
are in millimeters except volume, which is in cubic 
millimeters. 

4.2 A B-spline Model 

Given measured data points on the surfaces of the 
calibrated sphere and the two phantoms, surface mod- 
els can be constructed using B-splines as basis func- 
tions. In this section we will define cubic B-splines and 
show the construction of tensor products of B-splines. 



4.2.1 Cubic B-Splines in One Variable 

Suppose that a function;^ ^f(^), ^ ^ [^, b], is known 
at the n points (x^, yi) ,'", (x„, j^J, where a<Xi<X2< '" 
<x^< b, y^ =f(Xq), q = I,---, n. a and b are finite inter- 
val bounds. It is known that a polynomial of degree 
n- I, P(x), can be constructed to pass through these n 
points. In the case of highly accurate data points this 
polynomial can be constructed to interpolate these n 
points by, for example, Lagrange polynomial or 
Newton divided difference algorithms. But, for large n, 
it is also well known that polynomials of high degree 
can produce unwanted oscillations between the interpo- 
lated points. It is crucial then to approximate sets of 
data with as low degree polynomials as possible. Of 
course these polynomials may or may not interpolate 
the data points but may be made to come as close to 
them as possible. The ability to create highly flexible 
approximants from low-degree polynomials is a signif- 
icant advantage of functions called splines. 

Given a set of r real numbers, satisfying a< ^i< ^2 
< " < ^p<b, a spline function, F(x), of order p (or 
degree p- \) with knots ^i, ^2^'", ^p is di function that 
satisfles two properties: (1) in each of the intervals 
x<(gi; ^j_^<x<^j^ (7 = 2,3,-, r); ^p<x, F (x) is a 
polynomial of degree p- I or less. (2) F (x) and its 
derivatives of orders 1,2, •••,/> -2 are continuous. 
This would mean, for example, a spline function of 
order four would be constructed from polynomials of 
degree three (cubic) or less on the intervals x< ^^i ^j_i 
<x< ^j, (/ = 2,3,---, r); ^ p<x with continuous deriva- 
tives of orders one and two. 

A well known mathematical technique to construct 
complicated functions is to form linear combinations of 
simpler functions, called basis functions. In the current 
paper, data sets will be approximated by linear combi- 
nations of special spline functions, called B-splines, for 
Basis splines. It is known that any spline function can 
be represented in terms of B-splines. This particular 
basis has the advantage of leading to computational 
algorithms that are elegant, efficient, and stable. Only 
B-splines of order four will be considered here. They 
are cubic splines that are non-zero only over four 
adjacent intervals between knots (see Fig. 7). The nota- 
tion for a B-spline is Npj{x), where Npj{x), is 
zero everywhere except in the range 4_^<x<4, 
where in this work p = A. To simplify notation let 
Ni{x) = N^j{x). Then a B-spline of order four, or 
cubic B-spline, is a cubic spline with knots <^/_4, ^/-3, 
^i-2^ <^/-i5 ^i that is zero everywhere except in the 
range (^,_4<x<4. N^ix) is defined uniquely except 
for a scaling factor and is conventionally taken to be 
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&spline of Order 4 




Fig. 7. B-spline of Order 4 with knots at x = 0.1, 0.2, 0.3, 0.4, 0.5. 



positive throughout the range ^i_^<x< ^^ and has a 
single maximum value. Since N^ix) is a cubic spline 
it has continuous derivatives of order one and two at 
(^, _4 and 4 . These derivatives are zero at the endpoint 
knots. 

To define a complete set of B-splines on the set of 
points a< ^x< ^2^'" ^ ^r^b it is necessary to intro- 
duce eight additional points at the boundaries given by 

S-3? S-2? S-l5 So? Sr+l5 Sr + 2? Sr + 3? Sr + 4- ^^ IS USUal tO 

have ^Q= a, ^f^i = b. With this augmented set of knots 
one can define r + 4 fundamental cubic B-splines, 
Ni(x), / = 1, 2, •••, r + 4. Then the general cubic 
B-spline has a unique representation in the range 
a<x<b of the form 



These may be written in matrix notation as 

Ac^y, 



(19) 



F(x) = '^c,N,(x) 



(17) 



f=i 



where r=r+ 4. 

There are computational advantages in using cubic 
B-splines. For any given x, all but four adjacent Ni(x) 
are zero. In particular, if x g [(^^ _i, <^, ], the four non- 
zero cubic B-splines are Ni(x), Ni+i(x), Ni+2(^), Ni+3 i^)- 

A least squares curve fitting problem to the data 
a<Xi<X2< '" <Xj^< b, yq=f(Xg), q = l,"',n, becomes 
one of determining the coefficients c^ as the least 
squares solution to the equations 



r 



n 



(18) 



where A is the nxr matrix whose element in column / 
of row q is Ni(x^) and c, y are column vectors with 
elements c,, j^^, ^ = 1, •••, ^, respectively. If the data 
points are arranged in increasing order of x, then the 
matrix A becomes a banded matrix with bandwidth 
four. For a more thorough discussion of B-splines and 
their computation see de Boor [11] and Cox [12]. 

There are two functions in the MATLAB Spline 
Toolbox that implement the evaluation of B-splines. A 
set of knots can be augmented at the ends by the func- 
tion augknt and the B-splines and their derivatives can 
be evaluated by the function spcol. 

4.2.2 Tensor Products of Cubic B-Splines in Two 
Variables 

In this section the B-spline concept will be extended 
to two dimensions in order to fit two dimensional 
scattered data by a surface function. In this surface- 
fitting problem data points {x^.y^, q= 1,2, •••,^, and 
values at these points, z^ =/(x^, j^^), q= 1,2, •••,^, are 
given. The surface model used to fit these data points 
will involve sums of products of B-splines. 

To introduce this model, define a rectangle, say R, by 
a<x<b, c<y<d in the {x,y) plane. The definition 



/=i 
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does not restrict the data points to the EucHdean plane. 
They could, for example be angular coordinates, as will 
be seen later. The rectangle is subdivided by sets of 
knots ^1, i= 1,2, •••, r and rjj, J = 1,2,---, f , where 

a < (^1 < 4 < ••• < <^r < ^ and c < 7]i < 7]2 < ••• < 7]^ < d, 
where r, s are indices, not necessarily equal and 
a, b, c, d are the bounds on the rectangle R. The knots 
are then extended by eight in each dimension as done 
in the one dimensional case. These knots divide the 
rectangle R into rectangular panels in the plane given 
by Rij, i= 1,2, •••, r,j= 1,2, •••, s. Then, a basis set of 
splines for this pairing of knots can be constructed as 
products of B-splines Ni(x)Nj(y). In fact, the surface 
spline model is given by 



F(x,y) = yYq^N,ix)N.iy), 



i=l 7=1 



(20) 



called a tensor product of splines, where r = r + 4 and 

As in the one dimensional case, these tensor product 
splines have a number of advantages. First of all, the 
basis functions Ni(x)Nj(y) are each non-zero over a 
rectangle composed of sixteen panels in a 4x4 
arrangement. In particular Ni(x)Nj(y) is non-zero only 
when 4_^<x<4 and rij^<y<rij. Next, if (x^,j^^), 
q= 1,2, "',n and values at these points, z^ =/(x^, 3;^), 
q= 1,2, "',n diYQ given, then the fitting problem can be 
formulated as finding the least squares solution of 



7=1 ;=1 



^ = 1,2, 



n 



(21) 



i=l j=\ 



Again, this can be written in a matrix form as 

Ac ^ z , 



4.2.3 An Issue in Computing Tensor Product 

B-Splines and the Relation to Least Squares 

Unfortunately, for some choices of knots the result- 
ing matrix A in (22) might be rank deficient and it 
would not be possible to use the normal equations to 
solve the least squares problem. This can potentially 
happen in the case of widely scattered data, where some 
of the knot panels do not contain data points. This prob- 
lem can be solved, though, by using the matrix singular 
value decomposition. 

Assume that a tensor product spline has been com- 
puted as described in Sec. 4.2.2. A least squares fit can 
be done to the scattered data as follows. Since the 
matrix A could be rank deficient, with potentially zero 
rows or columns, we cannot rely on the standard nor- 
mal equations for determining the coefficients, but 
using the matrix singular value decomposition provides 
a convenient substitute (see Golub and Van Loan [15]). 
In the singular value decomposition of the matrix A, 
with the number of rows larger than the number of 
columns, the matrix is decomposed into the product of 
a column-orthogonal matrix U, a diagonal matrix S, 
with diagonal elements S^, and the transpose V of a 
square orthogonal matrix, so that A = USV. In order to 
solve the problem Ac = z ^q compute c = V(S^(U' z)), 
where 



s: = 




1 5^0 



5, 

otherwise 



(23) 



(22) 



defines the generalized inverse of S and U' is the trans- 
pose of U. A tolerance is used to determine which of the 
small singular values in the decomposition should be 
considered zero. 



where A is now a matrix with n rows and rs columns, 
and z is a column vector with n rows. The elements in 
row q, q= 1,2, •••, n, of A are formed as follows. We 
start with a fixedy,y =1,2, •••, s, and let /, / = 1,2, •••, r 
vary for thaty. Then, the element in column (/ - l)r + / 
for row q is given by N^ix^Njiy^ . They is then incre- 
mented and the /'s varied again. The end resulting 
column values in A for row q would look like 
Ni{Xq)N,{y^) N2{x^)N,{y^) - NXx,)N,{y^) N,{x^)N^{y^) - 
K(Xq)N2{yq)-N,{x;)N,{y;)-NXXq)N,{yq). For a full 
discussion of tensor products of B-splines see de Boor 
[11], Eberly [13], and Rogers and Adams [14]. 



4.3 A Knot Selection Algorithm 

In a least squares data fitting process involving 
B-spline basis functions, the resulting model residuals 
are sensitive to the knot placement for the B-splines. 
The selection of B-spline knots in order to achieve as 
small a residual as possible during the least squares 
process is a nontrivial task. One would be extremely 
lucky to manually select a set of knots that could 
achieve a very small least squares residual. In this 
section we will describe a heuristic algorithm that, in 
practice, generates a set of knots that produces small 
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least squares residuals. The strategy involves an itera- 
tive knot insertion algorithm. An initial set of knots is 
selected by the algorithm user and, at each iteration of 
the algorithm, new knots are inserted in the vicinity of 
the previous fit where the local residuals are the largest. 
The knots are not moved once they have been inserted. 
The iterative algorithm has a stopping criteria based on 
a statistical test. 

First we will discuss the knot insertion algorithm. It 
is based on a strategy suggested by Dierckx [3]. At the 
beginning of each iteration the assumption is that there 
exists a current set of knots. In the first iteration of the 
algorithm these would be an initial set chosen by the 
user. The tensor product of the B-splines is computed at 
the data points, a least squares fit is made to the data, 
and the current absolute residuals of the fit are comput- 
ed at each data point. The current knots divide the (x, y) 
plane into rectangles. Some rectangles have data points 
and others don't. Let the knots at the A:-th iteration be 
labeled a < (^/^^ < (^/^ < - < Q^^< b along the jc-axis 
and c < r\^^^^ < rif^ <-<r\,^^^<d along the j^-axis. The 
{k+ l)-iteration begins by associating all of the data 
points with the knot panels in which they fall. Let 
the index pair ij indicate the Rifih. panel defined by 
^i<x<^i+i, rij<y<rij+i. Then suppose the data 
values (jc/^\j^/^^), •••,(x/\j;/^), fall into the R,-th 
panel. Let Fj^ (x, y) be the least squares B-spline func- 
tion fit to the data at the A:-th iteration. We next compute 
the residuals of the fit at all of the data points 



Resid, =z -F,(x.v), q =1, 



ky q'> y q 



n 



(24) 



From these residuals we form the sums of squares of 
the residuals that fall within the knot intervals. The 
sums are separated into the x direction and the y direc- 
tion as follows. 

5^ =X{Resid^ : ^^ ^^^ <C}, i =U,-, r„ 

^•M{Residt:;7f <y,<C,}>7=U,...,.,. 

(25) 
What is meant here, for example in the case of 
5/t\ is that the sum, at the ^-th iteration of Resid J^, 

is taken over all pairs of points ( x^, y^ such that 
t^f <x^< t^lll . This is done for each / = 1, 2, • • •, ^ . 



We next find the maximums of these sums of 



squared residuals. 






1,2,- 
= 1,2, 



,''*}, 



■,■5* 



}> 



(26) 



where u = i for some / = 1,2,---, r^ and v =J for some 
J = 1,2,---, Sj,. The next step is to add one knot at a time 
at each iteration. In particular, a knot is added in the x 

direction, if 5*^^^ > S^!"^ or added in the y direction, 

' u.x v,y -^ ' 



if 5, ; < (5, , . The knot added is positioned based 

upon a weighted average of x or j^ data points in the 
columns or rows determined by the knot intervals with 
the maximum residual errors determined by (26). In 
particular, the positions are given by 



(^+1) 



ri 



(k+i) 



=1 



=1 



Resid, 

f^ „ . e(^) 



(k) 



u,x 



^q '' hu - \ ^ hu + 



(27) 



Resid 



kq 






yc 



: €' ^yq< Cl 



We note that 



I 



I 



Resid 



kq 



J(^) 

u,x 



Resid 



X : £f ' ^x< £*: \ = 1, 



1 



kq 



Xk) 



(k) 



v,y 



:li:'S>" <C', =1. 



U 



(28) 



SO that (27) represents weighted averages of all of the 
data pairs (x . v ) satisfying (^f ^< x < ^l^l\ and 4'^< v 



The knots are then reindexed as necessary and 
F^+i(x,j^) is computed as the least squares B-spline 
function fit to the data at the {k+ l)-iteration based on 
the new knot set. The iterations continue until the 
stopping criteria is met. 

The stopping criteria for the k-th. iteration used in this 
algorithm is based on the use of the 7^^-statistic, called 
the coefficient of multiple determination (see Draper 
and Smith [8]). R^ is the square of the correlation 
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between the vector of observed data, z^,q= l,"',n, and 
the least squares predicted data, Fj^Xq^y^, q= I,---, n. 
The statistic satisfies <R^ < 1 . This statistic is often 
used as a measure of how well the regression equation 
explains the variation in the data. It is known that in 
building models based on adding terms in the regres- 
sion equation, care must be taken in using this statistic. 
However, in the current algorithm, the statistic is used 
in a somewhat non-conventional manner. We use it as a 
measure of the benefit of adding more knots to the 
tensor product spline model. The knot selection 
algorithm is terminated when R^> 0.9^. 

Although the knot placement algorithm is heuristic, 
coupling it with the 7^^-statistic has shown good conver- 
gence in practice. It is reasonable to expect this, since 
knots are placed in intervals in which the fits at a pre- 
vious iteration showed the largest local error. Placing a 
knot there allows extra flexibility for those areas. 

4.4 Data Smoothing by Tikhonov Regularization 

As noted in Sec. 4.2.3, the data distribution can lead 
to a rank deficient least squares matrix. Even though a 
fit can be computed using SVD, evaluations of the 
fitted function at some points can lead to unwanted 
oscillations. It is necessary to introduce a smoothing 
procedure by regularization. Regularization is a way of 
introducing extra information to the least squares 
objective function in order to control overfitting of the 
data. It is the overfitting of data that can lead to the 
unwanted oscillations. In the present work a penalty 
term is introduced to control the smoothness of the 
resulting fitted function. The objective function will 
then balance the data fitting with the smoothness of the 
fit. The second partial derivatives of the tensor product 
B-splines will be introduced as the smoothing opera- 
tors. The objective function will take the form 



0(c) = £ 



q=l 



-l2 



yyc„N,(xjNJv)-z^ 






+1 

p=i 
t 

p=i 



n2 



r s 




Xc,.n;\u^)n.{v^) 



i=\ j=\ 



r s 



(29) 



-^2 




XcM,{u^)N;'{v^) 



i=\ 7=1 



evaluated at a new set of points chosen so that every 
knot panel has an equal number of points assigned to 
the panel. In the current case there will be a total of t 
points throughout the knot panels given by {Up,v^, 
p= 1, -,/. 

To write this in matrix form we will define two new 
matrices B^ and B2 as follows. For m= 1,2,"-, rs let 
ijmjm) be such that m = (j„- 1) r + /^ . Then define the 
matrix elements 



B,{p,m) = N,Ju^)Nl{v^). 



(30) 



We can rewrite (29) in the form 



rs 



n2 



0{c) = Y, Yi^,nAq,m)-z^ 

q=\ \_m=l 



rs 



n2 



p=\ {_ m=\ 



+ 






rs 



n2 



^X^-^2(A^) 



p-l L tn-l 

2 



(31) 



Ac-z 



+ 



XB,c 



+ 



XB^c 



U \ (z^ 



V 



XB^ 
XB. 



c — 



J 






V J 



The matrix dimensions in this objective function are: 



Ai^nv. rs, B^ is /^ x rs, B^is tx rs. 



V 



XB, 
XB. 



J 



is {n + 2t) X rs, and 



^z^ 








is(f2+2/)xl. 



If we perform a QR decomposition, then 



/ 



V 



A 

XB, 

XB. 



\ 



= QR. 



J 



(32) 



A is called a Tikhonov parameter. 

In the objective function the data values are given by 
{Xq,yq,z^, q=l,"',n. The smoothing terms will be 



where Qis(n + It) x(n + It) and Ris(n + It) x rs. Q is 
orthogonal, so that Q^ Q = I, and R is zero except in its 
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upper right comer. Let R be the rs x rs upper triangular 
portion ofR and let 



d = Q' 



^z^ 






V J 



(33) 



A 



Then, let d be the upper rs entries of d. The mini- 
mum of 0{c) then satisfies Rc = d. For a more 
complete discussion of the QR method see Golub and 
Van Loan [15]. 

At this point we need to discuss the selection of the 
Tikhonov parameter, X. A number of methods for 
choosing the parameter have been discussed in the 
literature (see [16, 17, 18, 19, 20]). For this work, 
though, we used a graphical method that led to the 
selection of a parameter in a few iterations. Before 
continuing to discuss the selection of the Tikhonov 
parameter for the current study we need to change the 
coordinates of the original probe points to spherical 
coordinates defined on a rectangle. 

4.4.1 A Coordinate Transformation 

Since our tensor product B-spline requires a rectan- 
gular grid, we convert our CMM data to spherical 
coordinates. We start with a given set of n points, 
(xi,yi,z^, i= 1,2, "',n, for some n, on a surface. As 
measured, these points are given relative to the CMM 
origin. The first step is to center the data by using the 
center-of-data mass of the data points. This establishes 
an origin at the center-of-data mass point. It is done in 
order to establish a common reference point interior to 
the measured data points. It also simplifies writing 
vectors from the origin to the data points and allows 
the introduction of spherical coordinates. The center- 
of-data point is given by 

- __Lv 

n ,=1 



1 " 
n i=i 

1 " 



(34) 



n 



i=\ 



Since the data set is enclosed in a near- spherical 
bounded region, it is reasonable to identify the 
Euclidean data points with spherical coordinates. This 
will map the points on the sphere to a rectangular 
surface where the surface coordinates are designated 
by 6, (j) and the height of a surface point is given by 
R((l), 0). In order to use spherical coordinates to repre- 
sent points on the boundary of a surface, we need to 
restrict our analysis to surfaces that are called 
star-shaped. These are surfaces in which a ray drawn 
from the center-of-data mass intersects the boundary 
in a unique point. Whereas, in the definition of the 
B-splines, we used the coordinates (x,y) we will 
now use (0 , 0) and build B-splines in terms of 
these spherical coordinates. Euclidean coordinates 
will now refer to the measured data points. This switch 
in notation should not, it is hoped, cause too much 
confusion. 

To each data point there is a vector from (x, y, z) to 
(Xi,yi,z;), given by F, = (x„ j^,, z,) - (x,>^, z). Further- 
more, any point within the bounding sphere can be 
identified by spherical coordinates of the form 
S(R,(I), 6) = (R sin(0cos(0), R sin(^)sin(0), R cos(^)), 
where x = R sin(^)cos(0), y = R sin(^)sin(0), z = 
R cos(O), for < ^ < 71, < < 2k. is referred to as 
the colatitude and (j) is referred to as the azimuth 
(Fig. 8). Table 2 associates the three dimensional 
octants with their spherical coordinates (where we have 
suppressed R). Therefore, the Euclidean coordinates were 
converted to 0, (j) angles on a rectangle [0, tt] x [0, 2n]. 
The height, ^(0, 6), at each 0, (j) is taken as the estimat- 
ed radius from the center-of-data mass of all of the 
Euclidean coordinates. 

As an illustration, the results of the conversion of the 
coarse probe points for the Green phantom from 
Euclidean coordinates to (j), coordinates are shown in 
Fig. 9. We note the density of data points near the equa- 
tor is higher than towards the poles at ^ = and =n. 
Unfortunately the distribution of data points was dictat- 
ed by the software controlling the CMM. The lack of 
data points near the poles leads to a well known prob- 
lem, called the Pole Problem in the literature (see 
Dierckx [3]). It creates rank deficient matrices during 
the least squares fitting process. Section 4.2.3 discussed 
how the singular value decomposition can be used to 
handle this problem. 
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Fig. 8. Spherical Coordinate angles. 6 is the colatitude and is the azimuth. 
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Fig. 9. Plot of the 6 , Coordinates of the Probe Points for the Coarse Data. 
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Table 2. Octant Equivalence between Euclidean and Spherical Coordinates 



Octant Number 



3-D Octants to Sperical Coordinates 

Cartesian Coordinate 



Spherical Coordinate 



1 

2 
3 
4 
5 
6 
7 
8 



x>0, y>0, z>0 

x<0, y>0, z>0 

x<0, y<0, z>0 

x>0, y<0, z>0 

x>0, y>0, z<0 

x<0, y>0, z<0 

x<0, y<0, z<0 

x>0, y<0, z<0 



0<e<n/2, 0<4><n/2 

0<e<K/2, Tll2<(p<Ti 

0<e<n/2, n <(t)<3n/2 
0<e<n/2, 3n/2<(t)<2K 
n/2<e<n, 0<(^<7r/2 
n/2<e<n, n/2<(t)<n 
n/2<e<n, n <(/)<37r/2 
n/2<e<n, 3n/2<(f)<2K 



4.4.2 Choosing a Tikhonov Parameter 

Once the original probe data had been converted to 
spherical coordinates, the iterative selection of the 
Tikhonov parameter for the current volume estimation 
problem proceeded as follows. First, a set of knots was 
selected with X=0 (i.e., no regularization terms) in 
order to produce an R^ > 0.98 as described in Sec. 4.3. 
We found that, for all of the data sets examined, only a 
very few extra knots were added beyond the initial 
set used to begin the knot selection process. This led 
to a rapid convergence to a set of fmal knots in 
all cases. These knots formed panels in the rectangle 
[0, 7i] X [0, 2k]. Next, nine points were uniformly 
chosen in each panel and the regularization terms 
formed. Numerical experimentation showed that the 
use of nine points in each panel provided sufficient 
extra data to the regularization terms in order to smooth 
the fmal fits. An initial Tikhonov parameter. A, was 
chosen and the objective function (29) was minimized 
by the QR method discussed above. An initial grid of 



40 6 and 80 (j) values was generated in the rectangle 
[0, Ti] X [0, 2k\. This grid was finally fixed on for 
selecting a Tikhonov parameter since further numerical 
experimentation showed that the response surface of 
the radii for denser grids did not change the final range 
of the radii. The coefficients, c, computed to minimize 
the objective function (29), were then inserted into the 
linear form Ac, where A was the tensor product matrix 
formed from all of the 40 x 80 grid points. The end 
resulting radii at the grid values were then computed 
and the spherical volumes for those radii were comput- 
ed. The maximum and minimum sphere volumes over 
the entire 40 x 80 grid were determined. To select the 
appropriate Tikhonov parameter, X values over a range 
were used and the differences between the maximum 
and minimum volumes were plotted. The value of A 
that produced the minimum difference was selected as 
the working A. For the current volume estimation prob- 
lem, the final A was A = 0.32. 

To illustrate the effect of using regularization terms 
to smooth the least squares fit of the radii see Fig. 10. 
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Fig. 10. Plot of Radius Data on a 40 by 80 Grid. A = 0. 
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Fig. 11. Plot of Radius Data on a 40 by 80 Grid. X = 0.32. 



This figure shows the radii data for the case of A = 0, 
or no regularization. Note the large radii oscillations at 
the boundaries (a range of approximately 10,000 mm). 
Now see Fig. 11 and note the narrow range (approxi- 
mately 0.1 mm) of the radii over the entire grid. This 
clearly shows the effect of the regularization terms on 
the least squares fit. The data used for these plots was 
the coarse Green phantom data. 



5. Volume Estimation by Surface 
Triangulation 

We could now estimate the volume using the triangu- 
lation method described in Sec. 2 by dividing the 
phantom surfaces into triangular surface patches as 
follows. First of all we partitioned the phantom 
surfaces at grid points located at the colatitude angles 
^1 = < ^2 ^ *"^ ^u ^ ^u+i ^ ^ from the north pole to 
the south pole and azimuthal angles 0i = < 02 < •*• < 
0;^ < 0;^ + 1 = 271 around the phantom surface, where v 
stands for vertical and h for horizontal. Since these grid 
points did not necessarily fall at the measured data 
points, a radius, 7^(0, ^), was calculated at the grid 
points from the regularized fitted B-spline surface 
model. At the north and south poles the radius was 
taken as the median value, Rnonh^ c)f the grid point 
values 7^(0^,0), / = I,---, h, for the north pole and ^south^ 
the median value of 7^(0^,71), / = \--,h, . The spherical 
coordinates of all of the grid points were converted to 
Euclidean coordinates on the surface by 



x = 7?sin(^)cos(0) ^<Q<n, 

y = R sin(^) sin(0) < < 2;r, 
z = RcosiO). 



(35) 



At this point we could apply the Divergence 
Theorem method of Sec. 2. The patches at the north 
poles were easily constructed to be triangular as part of 
the process of determining the contribution of each patch 
to the phantom volumes. In particular, at the north pole 
the designated point was (xi,j^i, Zi)= (0, 0, T^^^^^/^). 
We then iterated through the spherical coordinate points 
(0,, ^2)5 i^ l/'S ^- At each of these angle pairs there 
was a value ^(0^, 62), i= I,---, h. We generated the 
volume by adding up the contributions of each patch to 
the volume total. We did this by initializing a variable, 
vol, for the volume, to zero. We then started to generate 
the contribution from the first layer of patches at the 
north pole. As noted above, we set the north pole to 
(xi, j^i, Zi) and then set 

y^ = R((l)^,e^)sm(e^)sin((l)^), 
Z2 =R((l)^,0^) cos (6^), 



and 



X3 = R((l) 2,02) sin (62) cos((l) 2) 



2/' 



y, = R{(f> 2,62) sin {02)sm{(l) 2), 
Z3 = R((l) 2,92) cos (O2). 
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Sample Surface TrJangulation 



M 




-5 -5 



X 



Fig. 12. A Sample Triangulation of a Sphere Viewed from the North Pole (Vertex 1). 



These were set in a positive orientation. 

In order to make this process more concrete, we have 
included a sample surface triangulation with v= 5 
colatitude angles and h = ^ azimuthal angles in Fig. 12. 
In this figure the Euclidean grid points have been 
indexed. The North Pole (x^, y^, z^ is designated by the 
index 1. In the first step described above the points 
(x2,jF25-^2) ^iid (X3, j^3, Z3) are indexed in Fig. 12 by 
points 2 and 3 respectively. Next we computed the 
outward normal to the triangle patch by forming the 

vectors v^ = fe, j^2, ^2) - (-^b J^i, ^1), ^1 = (-^3. J'b. ^3)- 
(xi,j^i, Zi), for triangle 123 in Fig. 12 and then formed 
the normalized cross product 



n = 



v,xv^ 



(36) 



v,xv. 



We then computed the contribution that this patch made 
to the volume as 



V0l = V0l + (?2-(XpypZ^)) 
^ 3 



n 



V 



7=1 



i^jyj.i-k-^j.i-kyj) 



\ 



yj 



(37) 



where 



k = 



1<7<3 
3 7=3 



(38) 



is set in order to compensate for the cyclical vertex 
indexing around the triangle patch. Formula (37) is a 
combination of Eq. (6), where (n • (x^, y^, z^)) is the 
coefficient c in (6), and the second factor is a compact 
form of Eq. (11). The factors 1/3 from (6) and 1/2 from 
(11) were combined as a multiple of 1/6 after all of the 
summations had been performed. 
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We proceeded to the next patch in the North Pole 
layer. Again (x^, y^, z^ was the North Pole, indexed by 
1 in Fig. 12, and we then used the previous computa- 
tion to get X2 = X3, 3^2 ^ J^3? ^2 ^ ^3 ^^^ set 

X3 = R((l)^, 62) sin (62) cos ((j)^), 
Z3 - R((l)^, 62) cos (62). 

In Fig. 12 the new point (X3, y^, z^) is indexed by 4 in 
Fig. 12. The triangle of interest is now 134 in terms of 
indices. We computed the normalized cross product as 
for the first patch and then computed the contribution 
of the second patch to the volume using Eq. (37). We 
continued this process for 62, ^i, i= I,"', h. In Fig. 12 
we would have proceeded with computing contribu- 
tions to the volume by working through the indexed 
triangles 123, 134, 145, 156, 167, 178, 189, and 192. 

We next computed the contributions of the middle 
layer patches in a two step process. We iterated through 
each Oj, J = 2,---, v- \. The triangles at the South Pole 
were handled separately. For each Oj, J = 2,"-, v- \ 
and 0,, i=\,"',h, the patches were defined first in 
terms of four vertices to create four sided patches. Each 
of these patches was then divided into two triangles. 
The four vertices of a rectangular patch were identified 
counterclockwise as (Oj, 0,, R{i,j)), (^y+ 1, 0/, R{iJ + 1)), 
{Oj,,, 0,,i, R{i+\J+\)\ (Oj, (t),,,,R{i+\J)). As an 
example, in Fig. 12 one of the patches is identified by 
indices, in counterclockwise order, as 8 16 17 9. 
These indices would be associated with vertices 
(^2,01,^(1,2)), (^3,01,7^(1,3)), (^3,02,^(2,3)), 
(^2,02,^(2,2)). The four vertex patches were then 
divided into two triangles. For the first triangle in the 
rectangular patch we set 

Xi = R{hj) sin {Oj) cos (0,), 
yi = R{i,j) sin {Oj) cos {^^, 
zi= R(i, J) cos (Oj). 

X2 = R(i,J + 1) sin (^^+1) cos (0,), 
yi = R{Uj + 1) sin (^^ + 1) sin (0), 
Z2 = R{iJ+ l)cos(^^+i). 



x^ = R{i+ 1,7+ l)sin(^^ + i)cos(0 + i), 
3^3 = W + 1,7 + 1) sin (^^ + 1) sin (0.+ i), 
Z3 = 7^(/+ 1,7+1) cos (6/,.+ i). 



This triangle in the example Fig. 12 would be 
8 16 17. Again, the volume increment was computed as 
discussed previously. For the second triangle of the 
rectangular patch we maintained the same (x^, y^, z^ 
and set X2 = x^, y2 = y^, ^2 ^ ^3 ^nd then set 



X3 

yz = 

Z-i = 



R(i+ l,j) sin (Oj) cos ((^i+^X 

R(i+ l,y)sin(^^)sin(0.+i), 
R(i +hj) cos (Oj). 



This triangle in the example Fig. 12 would be 8 17 9. 
Again the volume increment was computed as before. 
We continued the process for Oj, j = 2,"-, v- \ and 
and 0, /= I,---, h. 

Finally, at the South Pole there were h triangles to 
include in the volume calculation. Their vertices were 
identified as follows 



Xi = R(i, v) sin (0^) cos (0), 
y^ =R(i,v) sin (0^) sin ((l)i), 
Zi =R(i,v) cos (Oy). 



X2 

yi 
^2 



0, 
0, 

-R 



south' 



X3 

y3 = 

Z-i = 



7?(/+l,i;)sin(^Jcos(0,+ i), 

7?(/+l,i;)sin(^Jsin(0,+ i), 
R(i+ l,v)cos(0^). 



As an example, in Fig. 13 one of these triangles 
would be indexed by the points 23 26 24, where the 
South Pole is point number 26, although the South Pole 
index may be hard to see in the figure. The contribution 
of these triangles to the volume was computed as 
above. After all of the triangle contributions to the vol- 
ume were computed, the final volume was taken to be 
vol = (1/6) vol. The factor 1/6 comes from the product 
of 1/3 from Sec. 2.1 and 1/2 from Sec. 2.2. The reader 
can refer back to these sections to see how the factors 
arose. The triangulation process can clearly be general- 
ized to a denser surface triangulation. 
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Sample Surface Triangulation 



r^ 




X 



5 5 



Fig. 13. A Sample Triangulation of a Sphere Viewed from the South Pole (Vertex 26). 



6. Computational Results 

Two computational processes were involved in 
generating the results for each object. The objects 
involved were the calibrated sphere and the two 
phantoms. First, sphere models were fit to an object. 
The second process involved three steps. The first step 
was to compute a good selection of knots for the tensor 
product spline function. The next step was to fit the 
tensor product spline function with the objective func- 
tion modified by the regularization terms to the object. 
Finally, once the tensor spline function had been devel- 
oped they were used to generate data at the vertices of 
the surface triangulation and to compute the volume 
using the Divergence Theorem. 

6.1 Computational Results for the Spherical Model 

Tables 3, 4, and 5 present the results obtained by 
fitting spherical models to the data from the calibrated 
sphere and the two phantoms. For the calibrated sphere 
there were five repeats for each of the coarse and dense 



data sets. Since the results of the point location repeats 
differed only at the micrometer level the data sets were 
combined by averaging to form two data sets represent- 
ing the coarse and dense data sets for the calibrated 
sphere. Similarly, the point locations of the three posi- 
tion data sets for both the coarse and dense data sets for 
the phantoms differed only at the micrometer level. 
Therefore, by averaging of the three position data sets 
for coarse and dense data sets, two working data sets 
were formed for the Green and Pink phantoms. 
Since the output of the Nelder-Mead algorithm 
was a final polytope surrounding the minimizing 
parameter with the polytope, the final reported para- 
meters were selected as the median value of the vertex 
values. These are the first four entries in the tables: 
Center x. Center y. Center z, and Radius. The units are 
millimeters. The fifth table entries are the spherical 
volumes based on the median radius value in cubic 
millimeters. The radius residuals were computed as the 
absolute value of 



Resid^. = yj(x,-cy+(y^-c^)^ + (z-cy- c, 



(39) 
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The sixth and seventh entries in the tables give the 
mean value and standard deviation of these residuals in 
millimeters. The eighth through the eleventh entries are 
the expanded uncertainties for the estimated Center x, 
Center y, Center z, and Radius, also in millimeters. 

The results in these sections on spherical model fit- 
ting involved the non-linear model (39) fitting to the 
sparse data sets generated by the CMM. This process 
did not involve the B-spline algorithm. Therefore any 
differences between coarse and dense results were most 
likely the consequence of the point selection in the 
metrology of the artifacts. 

6.1.1 Results for the Spherical Model Fit to the 
Calibrated Sphere With Coarse and 
Dense Data 

Table 3 gives the results of the fit of the sphere 
model to both the coarse (122 points) and the dense 
(181 points) data sets for the calibrated sphere. In both 
cases the average of the five repeat data sets was used 
for the fitting process. As can be seen, the two volume 
estimates differ by approximately 0.04 %. The radii 
estimates differ by about 0.01 %. There is a slightly 
larger difference when the spherical fit results are 
compared to the volume estimate based on the manu- 
facturer estimated calibrated sphere diameter of 
19.05 mm. This would lead to a volume estimate of 
3619.8 mm^ Since the method of estimating the dia- 
meter by the manufacturer was proprietary we could 
not independently verify the measurement. We could 
only compare it against our sphere model fit results. 
The manufacturer measured radius differs from the 
computed radii in Table 3 by about 0.08 % whereas the 
volume estimate based on the measured radius differs 
from the volume estimates in Table 3 by about 0.2 %. If 
we expand the radii in Table 3 to a diameter we find the 
values 19.0651 millimeters in the coarse case and 
19.0629 millimeters in the dense case. These differ by 
approximately 0.01 mm from the manufacture's meas- 
ured diameter. It would seem that these differences 
were sufficiently small so that the difference in value 
estimates for radius and volume were not considered 
significant. Although the estimated centers are slightly 
different, all other values are of the same order of 
magnitude. Based on these results we can accept that 
the CMM is producing accurate position data for spher- 
ical metallic artifacts. Tables 4 and 5, however, begin 
to show the consequences involved with attempting 
to measure slightly non- spherical and non-metallic 
artifacts with the CMM. 



6.1.2 Results for the Spherical Model Fit to the 

Phantoms With Coarse and Dense Data Sets 

From Tables 4 and 5 it is clear that the Green 
phantom is more spherical than the Pink phantom as 
indicated by the residuals and the expanded uncertain- 
ties. This simply confirms the fact that the Pink 
phantom was more difficult to measure using the vacu- 
um chuck due to its lack of sphericity. For example, the 
Mean Radial Residual for the Pink data is two orders of 
magnitude larger than for the Green data. The Standard 
Deviations for the Pink data are an order of magnitude 
greater, as are the Expanded Uncertainties. It is not 
clear how much the phantom material affected the 
results since it was difficult to set up a specific probe 
test for metallic versus phantom material. The artifacts 
would have to have been exactly the same size, 
positioned at exactly the right location, and probed at 
exactly the same coordinates to separate those factors 
from the material factor difference. There were no such 
comparable artifacts. We can make a guess, though, 
that there might be some effect due to probe force 
against the non-metallic material if we look at Table 3 
and Table 4 for the Green phantom, the most spherical 
of the two phantoms. The Expanded Uncertainties for 
the sphere fit to the calibrated metallic sphere and the 
Expanded Uncertainties for the sphere fit to the Green 
phantom data differ by one to three orders of magni- 
tude. Since the repeatability of the CMM measure- 
ments is at the 1 |im level, it is likely then that materi- 
al difference had some significant affect on the differ- 
ence in the uncertainties. It is a conjecture that this and 
the non-spherical shape of the Pink phantom account 
for a large part of the differences between the Pink 
phantom Expanded Uncertainties in Table 5, the Green 
phantom Expanded Uncertainties in Table 4, and the 
calibrated sphere Expanded Uncertainties in Table 3. 
For the Green phantom the radii estimates in Table 4 
between the fits of the coarse and dense data sets is 
approximately 0.3 % and the volumes differ by approx- 
imately 0.8 %. Whereas, for the Pink phantom the radii 
estimates in Table 5 between the fits of the coarse and 
dense data sets are 0.4 % and the volumes differ by 
approximately 1.3 %. In both Tables 4 and 5 the 
expanded uncertainties for the dense data sets are 
approximately an order of magnitude better than the 
results for the coarse data sets. This suggests that the 
radius and volume estimates in the case of the dense 
data sets for both the Green and Pink phantoms are 
probably the more realistic values, at least in terms of a 
spherical model fit. 
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Table 3. Results of a Spherical Fit to Coarse and Dense Data for Calibrated Sphere 



Properties 



Spherical Fit Results for Calibrated Sphere 

Coarse 



Dense 



Center x 

Center y 

Center z 

Radius 

Est. Volume 

Mean Rad. Residual 

Stand. Dev. Rad. Residual 

Expanded Uncert. Center x 

Expanded Uncert. Center y 

Expanded Uncert. Center z 

Expanded Uncert. Radius 



0.2803 X 10" 

0.4881 X 10" 

0.6642 X 10" 

9.5326 

3628.41 

0.3231x10" 

0.2156x10" 

0.2711x10" 

0.2711x10" 

0.4316x10" 

0.2188x10" 



-0.1331x10"^ 
-0.2034x10"^ 
-0.2564x10"^ 
9.5314 
3627.14 
0.8720x10"' 
0.22302 X 10" 
0.2357x10"' 
0.2357x10"' 
0.3757x10"' 
0.1909x10"' 



Table 4. Results of a Spherical Fit to Coarse Data for the Green Phantom 



Properties 



Spherical Fit Results to Green Phantom 

Coarse 



Dense 



Center x 

Center y 

Center z 

Radius 

Est. Volume 

Mean Rad. Residual 

Stand. Dev. Rad. Residual 

Expanded Uncert. Center x 

Expanded Uncert. Center y 

Expanded Uncert. Center z 

Expanded Uncert. Radius 



0.8216x10" 

0.1532x10" 

0.1832x10" 

10.1121 

4331.3 

1.4802x10" 

0.5526 X 10" 

0.1926x10" 

0.1926x10" 

0.2123x10" 

0.1146x10" 



0.9172x10" 

0.4275 X 10" 

0.1723 

10.0850 

4296.48 

0.8450x10" 

0.3093 X 10" 

0.4757x10" 

0.4753 X 10" 

0.7261 X 10" 

0.3632x10" 



Table 5. Results of a Spherical Fit to Coarse Data for the Pink Phantom 



Properties 



Spherical Fit Results to Pink Phantom 

Coarse 



Dense 



Center x 

Center y 

Center z 

Radius 

Est. Volume 

Mean Rad. Residual 

Stand. Dev. Rad. Residual 

Expanded Uncert. Center x 

Expanded Uncert. Center y 

Expanded Uncert. Center z 

Expanded Uncert. Radius 



0.7110x10"^ 

0.1329x10"^ 

0.1739x10"^ 

9.7332 

3862.09 

0.2315x10"^ 

0.2127 

0.2686 

0.2693 

0.2972 

0.1602 



0.2495 X 10" 

0.3272x10" 

0.3719 

9.7752 

3912.54 

0.3223 X 10" 

0.8935 X 10" 

0.3818x10" 

0.3862x10" 

0.5832x10" 

0.2880x10" 
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Table 1 shows that the diameter of the Green 
phantom is approximately 20 mm and that the diameter 
of the Pink phantom falls approximately between 
20.0 mm and 18.4 mm. This would imply that the 
sphere model for the Green phantom should have a 
volume of about 4188.8 mm^ With the estimate of 
0.3 mm uncertainty using the calipers, this would place 
the Green volume in the range of 4380.1 mm^ and 

4003.1 mm^. The estimated volumes in Table 4 are 
within this range for both the coarse and the dense data. 
The Pink phantom should have a volume between 

4445.2 mm^ and 3104.8 mml Table 5 for the Pink 
phantom shows that for the coarse and dense data the 
volume estimates fall from about 3862 mm^ to about 
3913 mm^. These volumes are within the expected 
range of the caliper measurements. 

6.2 Computational Results for the B-Spline Model 

In this section we will describe the method used to 
estimate volume uncertainties and discuss the results of 
using a B-spline surface model and the Divergence 
Theorem to estimate the phantom volumes. A nonpara- 
metric method to estimate the volume uncertainties was 
chosen since there did not exist any "ground truth" 
values for the Green and Pink FDA phantom volumes. 
All of the expanded uncertainty results are displayed in 
row two of Tables 6, 7, and 8. 

For each data set we conducted twenty one estimates 
of the volumes and volume uncertainties by increasing 
the density of the surface triangulation. The twenty one 
were selected since beyond that point the computer 
time became large and the volume increment per case 
was minimal. In order to plot the volume results in 
terms of grid density we use the term grid size by which 
we mean the product of the number of values times 
the number of (j) values for a particular grid density. As 
shown in Tables 6 through 8 we begin with a grid of 
ten and twenty (j) values and compute the volume and 
uncertainty. The grid size in this 1 x 20 grid is then 
200. The lvalues were then incremented by ten and the 
(j) values by twenty for each grid case until the last case 
of 210 values and 420 (j) values, giving a grid size of 
88200 (8.82 x 10"^). The volumes in Figs. 14 through 19 
grow very rapidly and appear to approach a fixed value 
as the grid size increases. They simply reflect the 
volume data in the Tables 6 through 8. 



6.2.1 Estimating Volume Uncertainties 

The uncertainties were estimated by the nonpara- 
metric "bootstrap" method. It is a computer intensive 
technique for estimating uncertainties and it involves 
repeated Monte Carlo resampling from the spherical 
coordinates of the original measured data sets with radii 
values modified by the fitting residuals, refitting the 
model to estimate new volumes, and finally computing 
an uncertainty for the process from the set of computed 
volumes. For a full discussion of the bootstrap see 
Efron and Tibshirani [21] but we will give a brief 
description here of how we applied the bootstrap 
method in the current study. 

The object of our application of the bootstrap was to 
develop volume uncertainty as a function of grid size. 
The process began with the conversion of the Euclidean 
data points to spherical coordinates. An automatic 
selection of knots was done. These steps were done 
once for a given data set. It was assumed that the 
modifications of the data made during the bootstrap 
would be small and not affect the knot selection. 
During the first pass of the bootstrap algorithm the radii 
of the spherical coordinates were fit by a tensor product 
of B-splines. The predicted values and residuals from 
this initial fit were called the master predicted values 
and master residuals respectively. The computed 
volume was then put in a list of volumes that was added 
to in subsequent passes of the algorithm. The algorithm 
then iterated two hundred times as follows. In the first 
iteration the master residuals were sampled randomly 
uniformly with replacement. The resampled residuals 
were then added to the master predicted radii and a new 
fit was performed. The new computed volume was 
added to the volume list and the master residuals 
sampled randomly uniformly with replacement for the 
next iteration. The process continued for all two 
hundred iterations. The standard deviation of the 
volumes in the volume list was then used as an estimate 
of the volume uncertainty and the average volume was 
taken as the reference volume for the chosen grid size. 
We found that in all cases the expanded uncertainties 
remained approximately constant and that is why only 
upper bounds for the expanded uncertainties are report- 
ed in the tables below. 
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Table 6. Estimated Volumes for Calibrated Sphere Data 



Calibrate Sphere Data. Volumes vs Grid Size 

3 3 

Expanded Volume Uncertainties: Coarse < 0.6 (mm ), Dense < 0.5 (mm ) 

3 3 

Coarse Volume (mm ) Dense Volume (mm ) 



Grid Case 








1 


10 


20 


2 


20 


40 


3 


30 


60 


4 


40 


80 


5 


50 


100 


6 


60 


120 


7 


70 


140 


8 


80 


160 


9 


90 


180 


10 


100 


200 


11 


110 


220 


12 


120 


240 


13 


130 


260 


14 


140 


280 


15 


150 


300 


16 


160 


320 


17 


170 


340 


18 


180 


360 


19 


190 


380 


20 


200 


400 


21 


210 


420 



3455.05 
3587.92 
3610.80 
3618.63 
3622.17 
3624.07 
3625.21 
3625.95 
3626.45 
3626.81 
3627.08 
3627.28 
3627.43 
3627.56 
3627.66 
3627.74 
3627.81 
3627.86 
3627.91 
3627.95 
3627.99 



3453.90 
3586.74 
3609.61 
3617.35 
3620.88 
3622.83 
3623.92 
3624.66 
3625.23 
3625.59 
3625.86 
3626.06 
3626.21 
3626.34 
3626.44 
3626.52 
3626.59 
3626.64 
3626.69 
3626.73 
3626.77 



Table 7. Estimated Volumes for Green Phantom Data 



Grid Case 


e 


(^ 


1 


10 


20 


2 


20 


40 


3 


30 


60 


4 


40 


80 


5 


50 


100 


6 


60 


120 


7 


70 


140 


8 


80 


160 


9 


90 


180 


10 


100 


200 


11 


110 


220 


12 


120 


240 


13 


130 


260 


14 


140 


280 


15 


150 


300 


16 


160 


320 


17 


170 


340 


18 


180 


360 


19 


190 


380 


20 


200 


400 


21 


210 


420 



Green Phantom Data. Volumes vs Grid Size 

3 3 

Expanded Volume Uncertainties: Coarse < 12 (mm ), Dense < 8 (mm ) 

Coarse Volume (mm ) 



Dense Volume (mm ) 



4125.34 
4283.92 
4310.93 
4320.67 
4324.90 
4327.17 
4328.53 
4329.41 
4330.01 
4330.44 
4330.75 
4330.99 
4331.07 
4331.33 
4331.45 
4331.54 
4331.62 
4331.69 
4331.75 
4331.80 
4331.84 



4092.77 
4250.85 
4277.65 
4285.66 
4289.86 
4292.12 
4293.47 
4294.35 
4294.95 
4295.37 
4295.69 
4295.92 
4296.11 
4296.26 
4296.38 
4296.47 
4296.55 
4296.62 
4296.68 
4296.72 
4296.77 



171 



Volume 115, Number 3, May- June 2010 

Journal of Research of the National Institute of Standards and Technology 



Table 8. Estimated Volumes for Pink Phantom Data 



Grid Case 


e 





1 


10 


20 


2 


20 


40 


3 


30 


60 


4 


40 


80 


5 


50 


100 


6 


60 


120 


7 


70 


140 


8 


80 


160 


9 


90 


180 


10 


100 


200 


11 


no 


220 


12 


120 


240 


13 


130 


260 


14 


140 


280 


15 


150 


300 


16 


160 


320 


17 


170 


340 


18 


180 


360 


19 


190 


380 


20 


200 


400 


21 


210 


420 



Pink Phantom Data. Volumes vs Grid Size 

3 3 

Expanded Volume Uncertainties: Coarse < 28 (mm ), Dense < 18 (mm ) 

Coarse Volume (mm ) 



Dense Volume (mm ) 



3668.27 
3811.59 
3836.22 
3844.59 
3848.41 
3850.46 
3851.69 
3852.48 
3853.03 
3853.41 
3853.70 
3853.91 
3854.08 
3854.22 
3854.32 
3854.41 
3854.48 
3854.55 
3854.60 
3854.64 
3854.68 



3719.71 
3863.30 
3887.99 
3896.39 
3900.21 
3902.27 
3903.51 
3904.30 
3904.85 
3905.23 
3905.52 
3905.74 
3905.91 
3906.04 
3906.15 
3906.24 
3906.31 
3906.37 
3906.42 
3906.47 
3906.50 



fO 



3840 
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Volume Estimates for Calibrated Sphere Coarse Data Set 



3 4 5 G 

Number of Grid Points 



8 



KlO 



Fig. 14. Volume Estimates for Coarse Calibrated Sphere Data Set vs. Grid Size. 
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Fig. 15. Volume Estimates for Dense Calibrated Sphere Data Set vs. Grid Size. 
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Fig. 16. Volume Estimates for Green Coarse Data Set vs. Grid Size. 
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Fig. 17. Volume Estimates for Green Dense Data Set vs. Grid Size. 
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Fig. 18. Volume Estimates for Pink Coarse Data Set vs. Grid Size. 
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Fig. 19. Volume Estimates for Pink Dense Data Set vs. Grid Size. 



6.2.2 Calibrated Sphere Volume Estimates and 
Uncertainties for the Coarse and Dense Data 

As a test of the accuracy of the B-spline volume 
estimation procedure, we applied the method to the 
coarse and dense data sets for the calibrated sphere. 
If the volumes computed by the B-spline method in 
Table 6 are compared with the sphere fit results in 
Table 3, we see that the B-spline method and the sphere 
fit method results differ by approximately 0.01 % for 
both the coarse and dense data sets. This indicates that 
the B-spline approach is performing as adequately as 
the straightforward sphere fit model and thus can be 
relied upon to produce good results when applied to the 
phantom data. 

6.2.3 Phantom Volume Estimates and 
Uncertainties for the Coarse Data 

In this section we discuss briefly the results of com- 
puting the volumes and volume uncertainties for the 
two phantoms using the coarse data sets. As described 
earlier, twenty one surface grid cases were used and the 
results are shown in Tables 7 and 8. Figures 16 and 18 
show the trends of the volume estimates as the grid 
sizes increase. The estimates rise rapidly, and both tend 
to what appears to be stabilized values. The figures 
graphically display the values in the tables. 



The spherical fit to the Green phantom data in 
Table 4 indicates an estimated volume of 4331.3 mm^ 
Table 7 indicates an estimated volume of 4331.84 mm^ 
for the largest grid size of 210 by 420 (j). This is 
approximately a 0.01 % difference suggesting that the 
Green phantom is very nearly spherical. 

The spherical fit to the Pink phantom data in Table 5 
indicates an estimated volume of 3862.6 mm^ Table 8 
indicates an estimated volume of 3854.68 mm^ for 
the largest grid size. The difference is approximately 
0.21 %, suggesting the slight non-spherical shape of the 
Pink phantom. 

From Tables 7 and 8, the uncertainties are shown 
bounded on the coarse surface grid by 12 mm^ for the 
Green phantom and 28 mm^ for the Pink phantom. This 
larger uncertainty for the Pink phantom might be due to 
the non-spherical nature of the Pink phantom. As indi- 
cated here, the uncertainties for the Pink phantom are 
approximately twice those of the Green phantom. 

6.2.4 Volume Estimates and Uncertainties for the 
Dense Data 

In this section we report the results of volume esti- 
mates and the extended uncertainties for the dense data 
sets for both phantom surfaces. The first thing of inter- 
est is that the stabilized volume estimates for the Green 
dense data in Table 7 differ from the volume estimate 
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for the Green coarse data by about 0.82 %. This is in 
the negative direction in that the volume estimates for 
the Green dense data were smaller than for the Green 
coarse data. In the case of the volume estimates for the 
Pink dense data in Table 8, the differences with the 
volume estimate for the Pink coarse data is approxi- 
mately 1.33 %. This is in the positive direction in the 
sense that the volume estimate for the dense data is 
slightly larger than the volume estimate for the coarse 
data. These results are consistent with the results of the 
spherical model fit to the Green and Pink phantoms in 
Tables 4 and 5. Another thing to notice is that the uncer- 
tainties are lower for the dense data sets as shown in 
Tables 7 through 8 compared to those estimated from 
the coarse data, although the uncertainty for the Pink 
phantom is about twice that for the Green phantom. 
These uncertainties are consistent with the uncertain- 
ties computed during the fits of the sphere models in 
Tables 4 and 5. It is not clear why the results for the 
dense data differ slightly in the directions they do from 
the results for the coarse data. A complete analysis of 
these differences is beyond the scope of this paper, but 
would be appropriate for a detailed study related to an 
analysis of the fitting algorithms and their sensitivity to 
the surface data distribution. 



7. Summary 

The B-spline surface model joined with the 
Divergence Theorem seems to be a viable approach for 
estimating volumes of the near-spherical molded 
phantoms, but the volume results seem to depend 
strongly on the distribution of the data points on the 
surface. In the case of sparse data the problem of 
extending a surface fitted to the data to a grid on the 
surface leads potentially to unwanted oscillations in the 
neighborhood of some of the sparse data points. In the 
current algorithm this problem was dealt with by 
adding regularization terms to the objective function. 
These provided a balance between the fit and the 
surface smoothness in order to control overfitting of the 
data. The final regularization parameter was A = 0.32 
for all data sets. 

The results obtained from both the coarse and dense 
data distributions appear to be consistent with the 
results from the spherical model fits. As shown in 
Table 7 the volume estimates for the Green phantom 
are very close to the estimates obtained by the spheri- 
cal model in Table 4 suggesting that the Green sphere 
is essentially spherical in shape. The volume estimates 
in Table 5 for the Pink sphere are higher by about 



0.19 % for the coarse data and 0.15 % for the dense 
data than the B-spline model estimates in Table 8 
suggesting that the Pink phantom might be close to 
spherical, although it does show some non-spherical 
tendencies. The uncertainties in the Pink phantom case 
suggest though that the Pink phantom may indeed have 
a slightly non-spherical shape. This is also suggested by 
the difficulty with holding the Pink phantom in the 
vacuum chuck during one of the surface data probing 
experiments. 

In summary we conclude that the combination 
of using the CMM to obtain surface data and the 
B-spline/Divergence Theorem volume estimation 
method can produce useful results but that there are 
some limitations. First, the CMM is primarily used for 
probing manufactured metallic artifacts and its use in 
probing non-metallic artifacts such as the FDA phan- 
toms can likely lead to some larger uncertainties in 
volume estimation. Second, the distribution of probed 
points can lead to results that indicate the models used 
are sensitive to the point locations. Third, the methods 
employed may not provide useful volume estimation 
for more complex artifacts that surely would arise in 
developing simulated lung cancer phantoms. The 
current artifacts were near- spherical and even these did 
not provide fully expected results based on data distri- 
bution for the Pink phantom. Further research on the 
affect of surface data distribution is required. Fourth, 
the Pole Problem is definitely a limitation that required 
regularizing of the objective function. Finally, the 
number of data points obtained on the surfaces was 
limited by the relative size of the artifact and the CMM 
probe size. 
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